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We present an analysis of the time evolution of self-interstitial atom and vacancy (point defect) 
populations in pure bcc metals under constant irradiation flux conditions. Mean-field rate equations 
are developed in parallel to a kinetic Monte Carlo (kMC) model. When only considering the elemen- 
tary processes of defect production, defect migration, recombination and absorption at sinks, the 
kMC model and rate equations are shown to be equivalent and the time evolution of the point defect 
populations is analyzed using simple scaling arguments. We show that the typically large mismatch 
of the rates of interstitial and vacancy migration in bcc metals can lead to a vacancy population 
that grows as the square root of time. The vacancy cluster size distribution under both irreversible 
and reversible attachment can be described by a simple exponential function. We also consider the 
effect of highly mobile interstitial clusters and apply the model with parameters appropriate for 
vanadium and a— iron. 

PACS numbers: 61.80.Az,61.82.Bg,61.72.Cc 



I. INTRODUCTION 

Many properties of metals depend crucially on the type 
and concentration of defects that perturb the ideal crys- 
tal structure. Of these, the simplest are point defects, 
such as self-interstitials and vacancies. Since their for- 
mation energies are of order electron volts, their equi- 
librium concentrations tend to be very low. They do 
form in abundance, however, in radiation environments 
due to the collisions between the irradiating species (elec- 
trons, heavy ions or neutrons) and the atoms of the host 
crystaii. When the energy of the impinging particle is 
close to, but above, the displacement threshold, the col- 
lision typically produces a single Frenkel pair. Particles 
with higher kinetic energy, for instance neutrons pro- 
duced in fusion reactions, create collision cascades which 
can produce not only Frenkel pairs, but an ensemble of 
mobile and immobile self-interstitial and vacancy clus- 
ters of different sizes. The ensuing evolution of the point 
defect distributions due to diffusion^ determines the long 
time scale degradation of the mechanical properties of the 
material, in addition to volume swelling at large doses. 

For body centered cubic (bcc) metals, such as a-Fe, V 
and Mo, molecular dynamics simulations have revealed 
a detailed microscopic picture of the diffusive motion of 
individual defects and vacancies'*. While the migration 
barrier AEy for vacancy migration is rather high (^ 0.5 
eV) , both self-interstitial atoms and small self-interstitial 
cluster are highly mobile and easily diffuse along partic- 
ular crystallographic directions ((lll)-directions). In V, 
the lowest energy self-interstitial configuration is a (111)- 
oriented dumbell, which migrates easily with a crowdion 
transition stated'. However, this easy migration, with 
barriers as low as AEi ~ 0.02 eV^, leads to long one- 
dimensional self-interstitial diffusional trajectories. The 
self-interstitial dumbbell can change direction by rotat- 
ing into other (lll)-dircctions. The barrier associated 
with such rotations AEr is of the same order as AE^ 
(i.e., AEr > AE,). 



In a-Fe, by contrast, the ground state of the dumb- 
bell self-interstitial is the (110) configuration, and access- 
ing the easy-glide (111) configuration requires overcom- 
ing a qualitatively similar rotational barrier that sepa- 
rates the (110) from the (111) configurations^. Although 
the migration barriers for the crowdion mechanism in a- 
Fe are also believed to be very small (< 0.04 eV), the 
observed effective migration barrier (including rotation 
into the (111) orientation and migration along the (111)- 
direction) are higher than in Similar arguments 

apply to Mo and Ta. Both diffusion mechanisms imply 
that interstitial transport in bcc metals takes place in 
the form of long one-dimensional trajectories with occa- 
sional directional changes that become more frequent as 
the temperature increases. In fee metals such as Cu, for 
instance, simulations showed that self-interstitial diffu- 
sion occurs through much more conventional, isotropic 
diffusion mechanisms^''. In most metals, however, large 
separations in time scale between self-interstitial and va- 
cancy motion exist. 

The intent of the present study is to illuminate the con- 
sequences of the intriguing microscopic diffusion mecha- 
nisms in bcc metals on the evolution of the point de- 
fect population using simple, but readily generalizable 
models. Our models shall minimally include local self- 
interstitial/vacancy mutual annihilation and absorption 
into (unsaturable) sinks which, in a real system, are in 
the forms of dislocations or grain boundaries. Together 
with the continuous production of point defects, these 
processes form the fundamental events in metals under 
irradiation conditions. 

On the continuum level, point defect dynamics can be 
described by a set of kinetic master or rate equations 
that treat the point defects as a continuous density whose 
temporal evolution is governed by various gain and loss 
processes. These equations do not take into account spa- 
tial correlations and are only analytically tractable in the 
simplest situations. Although we focus the discussion 
specifically on bcc metals, the continuum theory makes 
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no reference to the underlying crystal structure and could 
be readily applied to fee metals as well. In many cases, 
single interstitials and vacancies combine with defects of 
the same type to form stable clusters. The evolution of 
the cluster size distribution and other microscopic vari- 
ables is more conveniently studied in a particle based 
model that explicitly represents defects and their dif- 
fusion mechanisms (translation, rotation) on a lattice. 
The competition of events occuring with different rates 
can then be followed using a kinetic Monte Carlo (kMC) 
scheme. 

Both master equation and Monte Carlo approaches 
have frequently been employed to illuminate the physics 
of damage evolution in irradiated materials. Starting 
from the elementary processes mentioned above^, steady- 
state rate equations have been used to study the effects of 
preferential absorption of self-interstitial atoms at sinks 
(absorption bias)S'S, the formation of interstitial clus- 
ter during the cascade phase (production bias)Si2iiSiii, 
and the one-dimensional motion of these clusters in bee 
metals^iiSiii. kMC models were first used to evolve the 
primary damage state obtained from ns-long MD simu- 
lations to macroscopic time scales, but have increasingly 
been used to study the evolution of defect structure dur- 
ing continuous irradiation in copperiiSii^ii4, vanadiumpi^ 
and iron^^. 

Although many treatments have included a high level 
of atomistic detail from the start, we begin by analyzing 
the simplest situation that includes only production of 
Frenkel pairs, point defect recombination and absorption 
at sinks in Section |nj Even though this simplified case 
has been studied many times before!, we shall see that 
the inclusion of specific features of bee metals leads to 
a surprisingly rich evolutionary picture. This approach 
allows us to gradually increase the level of complexity 
of the model in a verifiable, controlled manner. We in- 
troduce first interactions between vacancies that lead to 
vacancy cluster formation in Section IIIII and then inter- 
actions between interstitials in Section Hvl In particular, 
we shall discuss how interstitial cluster mobility affects 
vacancy cluster formation. While the discussion in the 
text is applicable to bcc metals broadly, we apply our 
findings to the specific cases of V and a-Fe in Section Ivl 



II. POINT DEFECT DYNAMICS IN SIMPLE 
SITUATIONS 

A. Atomistic details and kMC model 

The physics of radiation damage evolution is governed 
by the production of defects due to irradiation and their 
subsequent elimination from the population through dif- 
fusional processes. Typical values for the defect pro- 
duction rate aF range between 10~^ dpa/s, in ion ir- 
radiation experimentsii, and 10~^° dpa/s, for neutron 
irradiation^"*, where F is the irradiation flux and a is the 
cross-section. A "displacement per atom" (dpa) refers to 



the production of one Frenkel defect pair per lattice site. 

Self-interstitial transport is composed of two parts, a 
one-dimensional diffusive motion along one of the 4 dis- 
tinct (111) directions and dumbbell rotations from one 
(111) direction to another. The temperature depen- 
dence of the (one-dimensional) diffusivity Di is usually 
described through the Arrhenius form in most cases ex- 
cept for vanadium, where the unusually low activation 
barrier can lead to more complicated behavior at higher 
T (see Ref.ll). For T < 600K, however, the interstitial 
diffusivity is well described by the Arrhenius expression 
Di/a? — fiyiex-pl—AEi/ksT], where AEi is the activa- 
tion energy barrier, Vi an attempt frequency and correla- 
tions are expressed through the correlation factor f {f = 
1 when the diffusion is uncorrelated) . The temperature 
dependence of the rotation rate 7^, by contrast, is always 
of the Arrhenius form, 7^ = t'r exp [— Ai^r/fesT], where 
AEr the characteristic rotation barrier and Vr is an at- 
tempt frequency. Values for both Vi and Vr range between 
lO^^s^^ and lO^^s"^. Finally, vacancies diffuse three- 
dimensionally with rate Dy/a^ = exp [— AS^/ZcsT], 
where AEy and Ur are the activation barriers and at- 
tempt frequencies, respectively. The defects can undergo 
two basic reactions: recombination once an interstitial 
and vacancy defect are within a certain "recombination 
volume", or absorption at sinks. Typical dislocation den- 
sities in metals are of the order 10*^ — 10*"* m^^, which 
translates into dislocation sink densities (per lattice site) 
of - 10"^ - 10-6. 

Since the relevant time scales (ns-sec) for the evolution 
of the point defect distributions far exceed those of molec- 
ular dynamics simulations, we employ a coarse grained 
description, in which we represent vacancies, dumbbell 
interstitial configurations and sinks as pointlike objects 
that occupy ideal lattice sites of a bcc lattice with lat- 
tice parameter a. In order to mimic constant irradiation 
conditions, new pairs of defects (self- interstitial and va- 
cancy) are introduced randomly onto defect-free lattice 
sites with rate aF. The microscopic diffusion process is 
replaced by instantaneous hops of point defects to vacant 
neighboring lattice sites with rates Di/a^ and Dy/a^, re- 
spectively. A self-interstitial is constrained to forward- 
backward hops along one of the four (111) directions, 
but can also rotate to another (111) direction with rate 
7r. Vacancies diffuse isotropically. A self- interstitial or 
vacancy recombines if it finds itself next to a vacancy or 
sclf-interstitial, respectively, or is absorbed if one of its 
8 neighboring sites contains a sink. At each time step, 
an event is chosen according to its probability and then 
executed. Time is advanced according to the usual con- 
tinuous time algorithm.**, where the time increment is 
chosen from an exponential distribution. 



B. Basic rate equations 

The elementary processes described in the kMC model 
can, alternatively, be described within a rate equation 



3 



formalism. Before the advent of large scale computer 
simulations, this method represented the only viable the- 
oretical approach for simulating long time radiation dam- 
age evolution. The rate equations can be solved by di- 
rect numerical integration or direct analysis in limiting 
cases. In the minimal model discussed above, the time 
evolution of the number densities of the interstitals rii (t) 
and vacancies ny{t) is given by the coupled nonlinear 
equation 
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~dt 



— (jF KyUJiyTli KiUlyiTly KigUligTii 

— GF KyUJiyTli KiUJyiTly KysLOysTly (l) 



Defect pairs are added to the population at a rate pro- 
portional to the particle flux F and a cross-section a. 
Loss can occur through a diffusing interstitial recombin- 
ing with a vacancy with rate and a diffusing vacancy 
recombining with an interstitial with rate uj^i. The re- 
combination rates are all proportional to the diffusion 
constant of the moving defect, but depend on the di- 
mensionality of the diffusion process. Ki and k„ are di- 
mensionless capture numbers that represent the spatial 
extent of the defects and their effective (possibly long 
ranged) interactions. Losses can also occur through ab- 
sorption at sinks with rates ujis and uJvs and correspond- 
ing capture numbers and for interstitials and va- 
cancies, respectively. An alternative representation of 
Eqs. Q is to define sink strengths fc^ via the relation 
k^D^ — Kx^xi where the subscript x refers to any of the 
combination of indices used above. 

The encounter rates of the defects are given by the 
number of distinct sites visited by a random walker 
per unit time. Since the mean squared displacement 
(i?^) = I'^N for a random walk with step length I, the 
number of sites visited is s = [{R"^) /P]^^'^ ~ iV^/^ in one 
dimension. By contrast, a detailed analysis of random 
walks on three-dimensional cubic lattices shows that a 
random walker visits 0{N) distinct sites after N hops. 
For a given density of target sites n, the typical collision 
time Tc is given by the condition Dtc/o? = 1/n (3D) 
and {DrJaY^'^ = l/« (ID), fr om which we deduce the 
encounter rates 



and 
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In a similar manner, the capture numbers k are also af- 
fected by the dimensionality of the random walk. Assum- 
ing ideal spherical defects of linear dimension r, k ~ r j a 
for a 3D random walk in a mean-field approximation^, 
but in general, capture numbers can also depend on 
spatial fluctuations and dose. For a ID random walk, 
the scaling with defect size becomes much strongesSi, 
At ~ (r/a)^. 

These expressions apply to strictly ID or 3D random 
walks. As discussed in the Introduction, we encounter 
an intermediate case in bcc metals, where rotations in- 
terrupt ID random walks and lead to 3D trajectories. 



The encounter rate of this random walk must, therefore, 
be larger than the purely ID case. For a rotation rate 
7r, the random walker performs on average D I^tO^ hops 
along a particular direction before rotating into a new 
direction. There are of these segments during time 
T. The mixed 1D/3D collision time thus follows from the 
condition {D / ^^o^'^Y^'^lrTc ~ ^/n, which implies 



^1D/3D = LOzD\JlrO?- ID 



(3) 



The mixed encounter rate uj\d /a^ scales like the 3D en- 
counter rate lozd, but is reduced by the square root of the 
ratio of the number of rotations to hops. This reaction 
rate has also been derived in ref. 21, This work views 
the kinetics in the intermediate 1D/3D regime as an en- 
hanced ID reaction rate, but the resulting expressions 
agree up to numerical prefactors. These authors also 
showed that the size dependence of the mixed capture 
number k '-^ (r/a)^, and ref. l22l provides an interpolation 
formula between the limiting cases using a continuum de- 
scription. While the encounter rates and reaction kinetics 
of random walkers decrease when the dimensionality of 
the random walk changes from three to one, the diffusiv- 
ity does not, since the mean squared displacement of a 
random walk of N steps of length / is (R(A^)'^) = 
independent of dimensionality. 

Note that these encounter rates are derived under the 
assumption of collisions with a stationary target. The 
case of several colliding ID random walkers becomes 
equivalent to a 3D random walk because, from the rest 
frame of a given walker, the other walkers appear to be 
executing a 3D random walk. This case would be rele- 
vant for describing the collisions of interstitials with each 
other, but not with the vacancies or sinks. 

Inserting these encounter rates into Eqs. assuming 
mixed ID/3D encounter for diffusing interstitials, yields 
rate eqations of the form 



drii 



= aF - riiTiy {Ky \f^Dil (i' - KiDy/a^) 

- KisUsniy/pDi/a^ (4) 
= aF - n.iny{K,y^fi3Di/a^ - KiDy/a^) 

- KysUsnyDy/a^, (5) 

where /3 = ^r<^ jDi is a dimensionless ratio that describes 
the relative frequency of dumbbell rotations and diffu- 
sional hops. 



C. Simple scaling analysis 

We now specialize these results to the common case 
of bcc metals, where "frO^ /Di <C 1 and Di/Dy ^ 1. 
Before performing the kMC simulations and solving the 
full rate equations numerically, it is instructive to analyse 
the limiting behaviors of this system^^. Inititally, there 
are no defects in the metal, and only the first term in the 
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rate equations is important. In this regime (regime I), 
defect densities increase linearly with time, 



aFt. 



(6) 



Once a sufhcient density of defects has been pro- 
duced such that the encounter times between defects 
becomes smaller than the time between creation events, 
loss through recombination (2"'' and 3'"'' terms) becomes 
important'^'*. As we shall see below, typical parameter 
ranges lead to a situation in which defect recombination 
becomes important before sink loss. The balance between 
creation and recombination makes a steady state possi- 
ble: drii/dt = duy/dt = 0. Ignoring the sink terms, one 
readily finds for this Regime II 



aFa^ 



-1/2 



p-1/2 



(7) 



which only depends on the dimensionless ratio F = 
(v73-D»+-Pf) ^ rpj^^ crossover from Regime I to Regime 
II occurs at time ti/n = Fati/n ~ F"^/^, i.e., ti/n ^ 
n^^ {Di / a'^)~^ . Note that this scaling regime is bound 
from above by the following condition: if riy becomes so 
large that the time between interstitial-vacancy encoun- 
ters is on average shorter than the time between two sub- 
sequent rotations (i.e., if F < /3~^), the scaling of the en- 
counter time becomes that of a ID random walk. In this 
case, the scaling of the steady state defect density with 
F changes from n-^ = n^^ ~ r~^/^ to nf = n^^ ^ F~^/'^. 
This regime is unlikely to be of experimental relevance. 

Eventually, loss through sinks becomes important and 
the steady state Regime II ends. Clearly, there exists a 
terminal steady state (see below) in which all loss terms 
in Eqs. balance the creation of defects. If Di ^ as 
in the bcc metal case, however, the system will initially 
loose mostly interstitials and very few vacancies. This 
breaks the symmetry between interstitials and vacancies 
and ^ Tii. The second steady state is, therefore, 
reached through a transient Regime III with distinct scal- 
ing. Subtracting the two rate equations (and neglecting 
vacancy loss at sinks), 

d{ny - Ui) dn^ r-Di 

T7 —rr = i^isns\/ P^n,, (8) 

at dt 0^ 

where Ug denotes the sink density. This equation allows 
us to eliminate in the rate equation for n^,, so that 



—— ~ at' — . 

dt KisTis dt 



(9) 



Integrating Eq. (O, we obtain (to leading order) power- 
law growth of 7iy with time. 



(10) 



In the case of extremely large sink densities Us, the en- 
counter rates ujis between interstitials and sinks would 
not be 3D as assumed above, but ID-like. In this case, 



ris would have to be replaced by nj., but again, such 
high densities are unrealistic. The crossover time tn/m 
between Regimes II and III can be obtained from the 
condition F ~ KisnsaFtu/m/ Ky, which implies tu/m 
{Ky/ KisUsY^^ / {Di/a?). The crossover occurs at constant 
time, independent of the defect densities. 

If there are only sinks for one species of defect and no 
sinks for the other. Regime III will simply continue. In 
all other final steady state will occur when all 

loss terms are important. Setting again ^ — = 
as in Ref. IT^ we obtain a condition for steady state. 



(11) 



where we have introduced a third dimensionless ratio a — 
Di/Dy. The interstitial and vacancy populations will, in 
general, be different. The steady state values are the 
roots of the quadratic equation 
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(12) 



Since typically <ti 1, the linear term can be neglected 
relative to the quadratic term, and we obtain the scaling 
of the vacancy density in Regime IV, 



F 



-1/2 



(13) 



The crossover time im/iv follows again from the con- 
dition KisUsaFtiiyiy/Ky ~ (F/av^)~\ i.e., ^m/iv 
{a\/]3Ky/Kisns)^^'^/{Di/a?) is again independent of Ui or 

Uy. 

The volume swelling is proportional to the excess va- 
cancy population in this regime, S = Uy — rii. Note that 
the imbalance arises here due to the very different diffu- 
sivities of the two types of defects. Upon termination of 
the irradiation, both the n;- and ni,-populations would 
relax exponentially to zero in this simple model. 



D. Numerical integration and kMC 

A full solution of the rate equations ^ is only possible 
numerically. In the following, we show a series of such 
numerical and kinetic Monte Carlo simulation results 
for a representative choice of parameters Di — lOOODt,, 
-fr = O.OlA/a^, so that a = 1000 and /? = 0.01 (see 
Section^for typical experimental regimes). In addition, 
we set the density of sinks to Ug — 10~^. Figure^shows 
the evolution of the average interstitial and vacancy den- 
sity (symbols) based on the above parameters obtained 
from the kMC simulation and numerical integration of 
the corresponding rate equations (solid lines) using the 
expression Eq. Q for the mixed encounter rates. The 
dimensionless parameter F was varied by changing the 
defect creation rate Fa. As predicted by the scaling anal- 
ysis, four distinct regimes appear with increasing time. 
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tDi/a^ 



and vacancy (0), ri^, densi 

- ^3 



10'',10^ and 10*' (F 



FIG. 1: Self-interstitial (■), 
ties as a function of time for F — 10^ , 10^ 
increases from the top to bottom of the figure). Time is mea- 
sured in units of the inverse interstitial hopping rate a? /Di. 
The solid lines show the result of direct numerical integration 
of Eqs. Q with = Ki = Kis = Kvs = 21, and the symbols 
correspond to the results of the kMC simulations in a periodic 
simulation box. The straight solid lines have slope 1 and 1/2. 



After first rising linearly with time (Regime I), rii and 
ri„ reach the steady state plateau of Regime II. Once 
loss through sinks becomes important (Regime III), n„ 
increases as t^^'^ while Ui decreases. Finally, all curves 
reach the steady-state Regime IV, where and are 
given by Eqs. lfTT|) and ifT^ . At the highest defect densi- 
ties, the kMC simulations where not carried out into the 
final steady state regimes, because the computational ef- 
fort becomes very large. 

Excellent agreement between continuum theory and 
kMC model for the 3 largest values of F is achieved 
by adjusting all capture numbers to a single numerical 
constant. In the simple situation examined here, kMC 
and rate theory are completely equivalent. Since we 
have used Wid/sd in the rate equations, this agreement 
also validates the scaling argument for the mixed 1D/3D 
encounter rate. For the two smallest F-values, we ob- 
serve increasing discrepancies between rate equations and 
kMC. As discussed above, we expect a transition from the 
mixed 1D/3D encounter rates coiv to ID dominated en- 
counter kinetics when the interstitial density becomes so 
high that they typically collide with vacancies or sinks 
before rotating (F < 100). This transition is not faith- 
fully reproduced by the rate theory which assumes only 
the limiting cases of either ID or mixed 1D/3D encounter 
rate scaling, but is properly captured by the kMC which 
includes explicit 1D/3D trajectories. 

Note also that with increasing F, Regime II begins to 
shrink and Regime I crosses over directly into Regime 
III. This happens because as the defect densities decrease 
below the sink density, point defect loss at sinks will be- 
come dominant before recombination plays a significant 
role. Since the crossover time ti/m — Kgns/ K^crF only 
depends on the point defect production rate, a further 
decrease of the production rate will eventually lead to a 



direct crossover from Regime I to Regime IV. A smaller 
sink density would push the appearance of Regimes 11 
and III to higher values of F. 



III. VACANCY CLUSTER FORMATION 
A. Irreversible aggregation 

To this point, neglected interactions between vacan- 
cies. Experimental vacancy-vacancy binding energies Eh 
can be of order the single vacancy migration energy. This 
suggests that it is not unreasonable to expect vacancies 
to form stable clusters or microvoids once they meet. A 
relatively simple approximation for this situation is to 
consider the limiting case of "irreversible aggregation", 
which neglects any dissociation of a vacancy from a clus- 
ter and thus sets Ef, = oo. The realistic situation, which 
includes finite binding energies, can be viewed as an in- 
termediate case between irreversible aggregation and the 
case of Eb — of the previous section. 

In the "irreversible aggregation" model, vacancies bind 
to form stable divacancies, leaving the population of free 
vacancies, while interstitials recombining with a diva- 
cancy recreate a single mobile vacancy. Stable cluster 
grow or shrink under the influence of arriving vacancies 
and interstitials. Denoting the number density of clus- 
ters of size m as nc{m), we can write the following set of 
rate equations^-: 



drii 



:(m) J 



Jt 



- 2KyUJyv)ny - {KyLJiy - KfuJis{2))n'i 

(Tfi — 1 m, \ 

<^iic(m-l) ~ l^y '-^vc{m))n'V 

+ {l^T'^^^ic{7n+l) - l^7'^ic{rn))ni, (14) 



where Uy = ndl) is now understood to refer to the 
free vacancy density (or cluster of size 1). The addi- 
tional terms in the equations for the evolution of rii and 
Uy account for interstitial/ vacancy-cluster encounters, 
vacancy /vacancy-cluster encounters, vacancy-cluster nu- 
cleation and divacancy decomposition. 0Jic(m) = 
nc{m)^/j3D,/a'^, LOydm) = nc{m)Dy/a^ and ujyy = 
UyDy/a? denote the corresponding encounter rates. The 
last expression in p4|l describes the evolution of stable, 
immobile vacancy clusters of size m > 1 due to the dif- 
fusive arrival of interstitials and vacancies. The hierar- 
chy of rate equations (|14|l is, in principle, amenable to 
a semi-analytical treatment via numerical integration, if 
all constants are specified. However, such a solution 
requires termination of the set of equations at a finite 
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cluster size rrimax, i-S-, one imposes a boundary condi- 
tion ncirriryiax) — 0. The corresponding kMC model does 
not require any ad hoc assumptions and takes all these 
processes naturally into account. 

Before examining this model in kMC, we can draw 
some premliminary conclusions. In Regime I, only nu- 
cleation of divacancies will be important. Since here the 
vacancy density grows linearly with time, we predict the 
scaling regime 



nc{2Y - {aFtf 



(15) 



The divacancy density grows cubically with time. In the 
steady-state Regime II, the vacancy cluster density is 
slaved to the free vacancy density and, therefore, is con- 
stant as well. For a ^ 1, we expect the total cluster 
density Uc = '"^cijn) to be much smaller than and 
rij, so that the presence of vacancy clustering does not 
yet strongly alter the interstitial and vacancy population. 
Using Eq. we find 



-1/2 



/(«V^)- 



(16) 



Once past Regime II, the cluster density will begin to 
rise again as n„ increases. n„ will then increasingly fall 
below the total number of vacancies in the system nytot, 
because increasing numbers of vacancies J2m>i iTT-ndm) 
become immobilized in clusters. These clusters act as 
additional sinks for the diffusing vacancies, but do not 
remove them from the system. The higher the cluster 
density, the smaller the loss of vacancies through sinks 
and interstitial recombination. We therefore expect that 
Regime III will extend to larger times and the presence 
of the immobile vacancy cluster delays the onset of the 
final steady-state Regime IV. 

For the total number of vacancies n^tot^ we can write 
the rate equations as 



dn„ 



dt 

dn- 



(17) 



m— 1 



m — l 



where we have neglected vacancy loss at sinks. This is 
the Regime III situation discussed earlier. Subtraction of 
Eqs. p7|l and elimination of Ui, as in Eq. Q, leads to 



driy 



dt 



aF 



K„ dn 



vtot 



^isT^s dt ^ ^ 
m—l 



-n^m). (18) 



At relatively early times, when only small clusters are 
present, I]m=i ^"c(to) — riytot and Uytot wiU initially 
increase as Uytot ~ {KisngaFty^'^ . The subsequent be- 
havior depends on the detailed form of the capture num- 
bers Kis- 

In the absence of vacancy clustering. Regime III ends 
when the vacancy density has become so large that va- 
cancy loss through sinks balances interstitial loss through 



sinks. This condition, Eq. Hll|l . must hold for a final 
steady state to appear. In the present simple model, 
this condition depends only on the ratio of the diffu- 
sivities, but in a real system, Eq. Hll|l also depends on 
sink concentrations and capture numbers, which may be 
different for interstitial and vacancies. Eventually this 
condition can also occur as a result of vacancy clus- 
tering, because vacancy clusters act as additional sinks 
for the interstitials and remove them symetrically from 
the system. At high cluster densities, the low mobile 
vacancy concentration can therefore be compensated. 
In the steady-state regime IV, we expect the hierarchy 

Because of the later en- 



< n. 



IV 



< n: 



IV 



< n. 



IV 



try into Regime IV, S = Uytot — n-i (the total amount 
of volume swelling) will be larger than in the absence of 
clustering. 

In Regime IV, we can obtain an approximate scaling 
relation for the total cluster density. The density is deter- 
mined by the competition between nucleation of clusters 
out of the vacancy gas and decomposition of divacan- 
cies due to arriving interstitials. Denoting the fraction of 
clusters that represent divacancies as /, we can write the 
rate equation 



dnc 



IV 



dt 



= KyulDy/a^ - f{t)i^^nfn,y/pD^/a', (19) 



which implies nj,^ — Kyn^/ /nfay^r, 



Kyny/fK2, 



where we have used Eq. Hll() . The cluster density is there- 
fore of order the free vacancy density with a numerical 
prefactor that depends on the mean cluster size. 



1. Average cluster density 

The effect of irreversible vacancy aggregation on the 
defect dynamics as seen through kMC is shown in Fig. [21 
for three values of F = 10^ F 10^ and F = 10^ with 
all other parameters as before. In this figure, we also 
plot the total cluster density Uc and the total vacancy 
density Uytot- In Regimes I and II, Uc and are nearly 
unchanged. The defect density Uc rises first ~ t^ and then 
remains constant during Regime II as discussed above. In 
Regime III, Uy begins to fall below the result from Fig.Q] 
because of trapping of vacancies into clusters. At the 
same time, Uytot rises with an ideal t^^^ law, and Regime 
III extends to larger times. Uc also rises again due to 
new cluster nucleation events out of the vacancy gas. All 
quantities reach constant values in regime IV. As one 
might expect, the total number of vacancies in the system 
is higher than in the absence of clustering. However, the 
general sequence of regimes remains unchanged. 

In order to gain additional insight into the growth dy- 
namics, we also numerically integrate the full set of rate 
Eqs. (solid hues in Fig. This requires specifica- 
tion of the capture numbers k™. In the simplest mean 
field model for nucleation dynamics, the capture numbers 
do not depend on cluster size m and „ = k™^ = const. 
(i.e., the point cluster model). Interestingly, we find 
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tDi/a^ 

FIG. 2: Interstitial densities rii ( ■), free vacancy density n„ 
(0) as well as the total vacancy density n^tot {*) and vacany 
cluster density ric (a) as a function of time for (a) F = 10^, 
(b) F = 10'', and (c) F = 10^. The thin solid lines show 
the result of direct numerical integration of Eqs. 1141 with 
~ = 25 and a rrimax ~ 20 (see text). 



surprisingly good agreement between the rate equations 
(solid lines in Fig. |2Jl and the full kMC simulations 
when using this simple approximation. This suggests 
that for relatively small m, is only very weakly size- 
dependent. The rate equations faithfully reproduce the 
sequence of regimes, but begin to show deviations at late 
times when larger clusters become more prominent and 
the point cluster assumption is no longer accurate. Here, 
begins to show some size dependence. The success of 
this comparison shows, however, that a mapping between 
kMC and rate equations is possible even with void nuclc- 
ation and growth if the k^'s are accurately parametrized. 
Note that in the present model, the volume swelling 



rate dS/dt is zero in Regime IV. This is possible be- 
cause we have chosen the same capture radii Kg for the 
sinks for both interstitials and vacancies, i.e., the sinks 
are symmetric with respect to defect type. S' ~ t would 
be valid, if, for example, the capture radius for intersti- 
tials is larger than that for vacancies. This situation can 
arise with the introduction of dislocation sinks (known 
as "dislocation bias" ) and has been invoked to explain un- 
usually large swelling rates (see also Section lVl|) at times 
later than considered hereSsS. 



2. Cluster size distribution 

A full description of the evolution of the point defects 
in the material includes a characterization of the clus- 
ter size distribution. The growth of these clusters is the 
result of the diffusive arrival of interstitial and vacancy 
defects at already nucleated clusters. The full dynam- 
ics of this process is described by the last equation in 
()14|l. Consider this equation in the following simplified 
notation: 



dnc{m) 
Jt 



n„(K™ ^nc{m - 1) - K'^nc(m))Dy /a'^ 
n,{K^+^n,{m + 1) - KTn,{m))^D,l i^Q) 



If rii and n„ change very slowly in time, the distribution 
nc{m) is given by the steady state solution to Eq. H2U|I . 
i.e., dnc{m) / dt — 0. This equation is supplemented by 
the detailed balance condition 

n,Kl^n,{m)Dja' = n^KT^^n,{m + l)^D,/a\ (21) 

which should also hold for m > 2. From this condi- 
tion, we can deduce the distribution ndm) by induction. 
Since riv = nc(l), we have nc(2) = n^Ki / {Kfa^/Prii) and 
consequently for all m > 1 



ndm) 



n 



1=1 



(av]3no™"' n: 



(22) 



1=2 ' 



Equation H22|l is tested against the kMC results in 
Fig. |31 where we show the cluster size distributions at 
four different times. All curves fall along straight lines in 
a semilogarithmic plot, and the slope decreases with in- 
creasing mean cluster size. A comparison with Eq. (|22|) 
requires, again, knowledge of the capture numbers k™. 
As in the previous section, use of the point cluster ap- 
proximation kI^/kI^ = 1 yields excellent agreement be- 
tween between the kMC data and Eq. H22|) upon inserting 
the values for ni{t) and ny{t) at the appropriate times. 

Note that the distributions shown in Fig. O are not 
peaked, i.e., the frequency with which different clus- 
ters appear decrease with increasing cluster size and sin- 
gle vacancies occur most frequently. This situation is 
in sharp contrast to other cluster growth situations as, 
e.g., found in submonolayer island growth during vapor 
deposition^^, where the distribution ndm) peaks at the 
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5 10 15 20 



m 

FIG. 3: Plot of the cluster size distribution nc{m) found 
from the kMC simulations with F — 10'' of Fig. HJa) at four 
different times tDi/a^ = 3 x 10** (■), 3 x 10^ (♦), 3 x lO'' 
(★), and 3 x 10^ (a) in Regimes III and IV. The straight lines 
correspond to n^t)"" /{a^n^{t))"^-\ 

mean cluster size (m). In this problem, however, clus- 
ters cannot shrink, since vacancies are absent. One ex- 
pects Eq. H22|l to hold as long as the assumption of quasi- 
stationary values for and riy is valid. 

B. Reversible attachment 




2 4 6 8 10 12 14 



The above discussion immediately raises the question 
of whether the results for ndt) and its size distribution 
Eq. (|22|l survive in the more realistic situation of re- 
versible cluster growth, i.e., vacancies have a finite prob- 
ability to attach to and to leave the cluster. Let us 
introduce "detachment rates" from a cluster of size m, 
Ideti'm) / for all vacancies, independent of their local 
environment. Equation H2U|I needs to be generalized by 
the addition of two terms, 

= n,(C"'«c(m-l)-K>,(m))i?,/a2 

+ n,(C+'«c(m + 1) - 

+ {n^+^n,(m + 1) - K^n,{m))^det{m) / a^A) 

where represent "detachment numbers" in analogy to 
the capture numbers k™,. Consequently, the condition 
for detailed balance now reads 

+ KJ';i+^n,{m,+ \Yidet{m)la^{2h) 
and Eq. (|^ generahzes to 

n^[rn) = ^^'-i ^""^ ^ (26) 

where a'{m) = jdet{'m)a'^ /Dy. We see that the general 
form of the distribution is the same, but the prefactor 
changes due to the additional growth and shrinkage prob- 
abilities. From a statistical point of view, interstitial ar- 
rival and vacancy detachment are equivalent. Eq. H2t)|) 



FIG. 4: (a) Interstitial density rii (■), free vacancy density 
riv ( 0) as well as the total vacancy density ((a)) and vacancy 
cluster density ric (*) as a function of time for F = 10^ , and 
10^ for reversible aggregation, where S — 0.01. (b) Cluster size 
distribution ndm) for the case of F = 10'' at four different 
times in Regimes III and IV. The straight lines correspond to 
n„(t)"7(aV^ni(f) + t^la')"'-'^ with = 0.25. 

can be easily evaluated for any functional form of the size- 
depedent capture numbers and detachment rates. In the 
point cluster approximation, where all >^Tvdl^\vd — 1 
and a' = a'{m) is size-independent, it predicts an expo- 
nential distribution as in the case of irreversible attach- 
ment, 

«cM = 7 7= " , r. (27) 

Figure 01 presents a numerical investigation of re- 
versible vacancy cluster growth using kMC. Although de- 
tachment rates may be cluster size dependent in general, 
we only introduce one size-independent rate ^detjo^ — 
^XyVD^ja^ for vacancy detachment for simplicity. This 
model would be most relevant for faceted cluster shapes 
with one dominant detachment rate from the faces. The 
other parameters are that of Fig. [3 We see in Fig. EJa) 
that the cluster density behaves in a qualitatively similar 
manner as in the irreversible case, but ric is reduced and 
the final steady state is reached at earlier times. As in 
the irreversible case, the cluster size distribution nj^m) 
shown in Fig. 01b) is well described by Eq. H27|) . The 
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FIG. 5: Interstitial densities ( ■), free vacancy density n„ 
(0), total vacancy density n^tot (a), vacany cluster density 
He {*) and density of immobile interstitial cluster (□) as a 
function of time for F = 10®. 



present discussion is of course only relevant when void 
coarsening can be neglected. 



IV. INTERSTITIAL INTERACTIONS 

One of the fascinating aspects of point defect dynamics 
in bcc metals is the high mobility of interstitial clusteraSi. 
The interaction between single interstitial atoms is even 
stronger than that between vacancies, and the intersti- 
tial clusters, such as dislocation loops^'*, are stable at 
all relevant temperatures. Unlike the immobile vacancy 
cluster, however, the interstitial cluster migrate easily for 
clusters with up to 50 or 100 interstitial atomsiSiSS. We 
now modify the preceding analysis to account for this 
behavior. 

In this analysis, we return to the irreversible aggre- 
gation limit, i.e., interstitials never separate after en- 
counter. This implies a reduction of the density of mobile 
interstitials due to nucleation of interstitial clusters, i.e., 
interstitials act as sinks for other interstitials. However, 
this creates a population of larger clusters with a larger 
crosssection and eventually larger capture numbers. The 
reaction kinetics will be affected by the competition of 
these two effects. 

The rate equations H14|l can now be expanded to in- 
clude terms representing the above processes, but be- 
come even more complex. In our kMC model, we begin 
by studying the effect of interstitial cluster formation by 
considering completely immobile interstitial clusters in 
complete analogy to the vacancy cluster. This situation 
is rarely realistic, but provides an upper bound on the 
magnitude of the effects. Figure |31 shows the evolution 
of the free interstitial and vacancy densities as well as 
the interstitial and vacancy cluster densities for the same 
parameters as before. As expected, the immoblization 
of interstitials increases the total number of vacancies 
in the system. In the final steady state, riytot is about 
three times as large as in the situation in Fig. |21 The 
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FIG. 6; Interstitial, vacancy, interstitial cluster and vacancy 
cluster densities for F = 10^ and F = 10^ (symbols as in 
Fig. 01. Interstitial clusters (□) diffuse and rotate with the 
same rate as single interstitials. The interstitial cluster den- 
sity first rises due to nucleation events, but then rapidly drops 
as mobile cluster collide with sinks. 



interstitial clusters (open squares) nucleate earlier than 
the vacancy cluster, but their density later drops below 
that of the vacancy clusters, since the single interstitial 
density is much lower than the single vacancy density. 

The situation in Fig. can be favorably contrasted 
with that of completely mobile interstitial clusters, i.e., 
the cluster diffuse with the same rates as the single in- 
terstitials regardless of their size. This case is again not 
fully realistic, since interstitial clusters seldom rotate into 
other (111) directions once they contain several inter- 
stitials (i.e., they diffusive one-dimensionally^^''^'^). For 
small clusters, however, the completely mobile interstitial 
cluster case represents a good approximation. Figure |H| 
shows the corresponding results for the various desities 
introduced above. As in the case of immobile intersti- 
tial clusters, the mobile interstitial cluster density first 
rises due to nucleation events. However, the highly mo- 
bile interstitial clusters also collide with the sinks, and 
the cluster density decreases rapidly. Since the free in- 
terstitial density also declines, nucleation events become 
rare. Once steady state has been achieved, the inter- 
stitial density has become so low that the nucleation of 
new interstitial clusters due to diffusion is almost com- 
pletely absent. Consequently, the vacancy densities and 
the total swelling rates are the same as for the case of 
non-interacting interstitials. 



V. APPLICATION TO VANADIUM AND IRON 

In this section, we apply our model to the V and a— Fe 
cases. Vanadium is particularly interesting, because here 
the effect of mixed 1D/3D diffusion is most pronounced. 
First principles calculations and classical MD simulations 
of V have yielded estimates of AEi = 0.018eV^ AEr 
0.44eV^ and AE^ w O.SeV*. For Fe, the different ground 
state of the self-interstitial ((110) instead of (111)) leads 
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to a higher effective activation barrier for ID migration, 
AE„ which ranges between 0.12 eV^" and 0.17 eV21. The 
rotation barrier was estimated as AE",. = 0.16 eV^, and 
the vacancy migration energy is assumed to be of the 
same order^^ as in V (the prefactors for all processes 
tend to vary by less than an order of magnitude). One 
therefore expects that the self-interstitial trajectories in 
a—Fe will be much more isotropic than in V. 

The two metals are best compared in terms of the rel- 
evant dimensionless parameter a, /3 and T. Table^sum- 
marizes the values for these quantities using the above 
energy barriers and the production rate for ion irradia- 
tion aF — 10~'^s~^ for two representative temperatures 
T = 300K and T = 600K. We first note that T is typically 
larger than 10^" and would in fact reach 10^° when typ- 
ical production rates for neutron irradiation are used. A 
kMC simulation with such large values of F would require 
very long simulation times, since diffusion and recombi- 
nation occur much more often than the introduction of 
new defects. It also requires large system sizes, because 
the defect densities become very small. In addition, the 
very small value of f3 implies very long ID segments of 
the interstitial trajectories between rotations. The pe- 
riod of the kMC simulation box should be several times 
larger than the typical length of those segments in order 
to properly reproduce the continuum theory values of the 
encounter rates. If the box period is much smaller than 
the ID segment, the trajectory will wrap around the box 
many times, but is not necessary space-filling. 

The last of these issues can be addressed by using a 
result from Section^ where we showed that the mixed 
1D/3D encounter rates scale like the ideal 3D rates that 
are reduced by a factor ^/p. We can therefore replace the 
explicit 1D/3D trajectories of the interstitials and inter- 
stitial cluster with ideal 3D random walks, but reduce the 
hopping rate by a factor of \/]3. This procedure leaves 
the reaction kinetics invariant (up to small corrections 
from the capture numbers) and implies that the effec- 
tive ratio of time scales between interstitial and vacancy 
migration is given by a' = Di^/P/Dy ~ yJUcwJa^ ■ In- 
terestingly, this ratio is very similar for both V and Fe 
at 300K and 600K (see Tabled, even though the values 
of a and /3 are very different. At a given temperature, 
the self-interstitial trajectories in a-Fe are much more 
isotropic than in V, but the encounter rates with vacan- 
cies and sinks only depend on the product of diffusivity 



TABLE L Dimensionless parameters a = Di/D-u, P ~ '^r/Di, 
V = {y/pDi + D^)/aF, and a' = for V and Fe at 300 K 
and 600 K. 
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FIG. 7: Interstitial, vacancy, interstitial cluster and vacancy 
cluster densities (symbols as in Fig.|SJl with parameters a' = 
10^ Us = 10-^ for (a) F = 10^ (b) F = lO" and (c) F = 
10*^ representative for V or Q-Fe at 600K. Thin solid lines 
show the result of numerically integrating Eqs. 1141 again for 
Km = 25 and rrimax = 20. The solid lines have slope 1 and 
1/2, respectively. 



Di and rotation rate 7^. Within the present model, one 
therefore expects that the point defect reaction kinetics 
in these two metals is very similar. 

In Fig. we use a ratio of time scales representative 
for V/a-Fe at T = 600K and show resuhs for F = 10^ 
F = 10^^ and F = 10^^. Here, we have multiplied the 
time axis with the production rate, which is a more con- 
vential presentation of the data in experimental stud- 
ies. Rescaled by the production rate, all curves initially 
coincide and start out with near-constant slopes. For 
the present parameters. Regime II is practically absent, 
and the vacancy density crosses over immediately from 
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Regime I (where it rises linearly with dose) into Regime 
III. The crossover dose is roughly constant and depends 
on the sink density. Note that the growth of riytot in 
Regime III is well described by a {aFty^'^ power law 
over several decades before the final steady-state Regime 
IV is reached. Since the onset of Regime IV is indepen- 
dent of time (see Sectional, larger values of F push the 
beginning of steady state to smaller doses. Here, steady 
state is reached at about 1 dpa for F = 10^ and earlier 
for F = 10^^ and F = 10^'^. There is some initial nucle- 
ation of interstitial clusters at F — 10^, but as discussed 
before, all of these rapidly disappear due to fast collisions 
with sinks. For higher values of F, the interstitial cluster 
density becomes negligibly small. 

In the same figure, we also show the predictions of a 
numerical integration of the full set of rate equations H14|l 
using the point cluster model (size independent capture 
numbers). As in Fig. [3 the agreement with the "exact" 
kMC simulation is satisfactory. Since the computational 
effort for direct kMC simulations for F > 10^° becomes 
tremendous, the rate equation approach, when properly 
parametrized, is clearly preferable. 



VI. CONCLUSIONS 

We have studied the early-time evolution of point 
defect populations with specific reference to migration 
mechanisms in bcc metals under constant irradiation con- 
ditions. Only very simplified models that incorporate the 
most important processes were employed in order to iden- 
tify the rate-limiting events. These models were solved 
using scaling arguments, direct numerical integration of 
kinetic rate equations and full kMC simulations. 

In the simplest case, which only considers homoge- 
neous defect production, recombination and sink absorp- 
tion, but no interactions between defects of the same 
type, the kMC model and the corresponding rate equa- 
tions were shown to be in near perfect agreement. We em- 
ployed simple random walk arguments to derive a mixed 
encounter rate that describes one-dimensional diffusion 
with occasional rotations. This encounter rate scales lin- 
early with target density as the isotropic 3D encounter 
rate, but is reduced by the square root of the ratio of 
rotation rate and hopping rate. As in ref. four dis- 
tinct scaling regimes for the point defect density with 
time where identified. First, the point defect densities 
increase linearly in time due to production (Regime I), 
then saturate as defect recombination sets in (Regime 
II). Sink absorption then begins to reduce the intersti- 
tial density and the vacancy density grows in time with a 
characteristic ^ t^^'^ behavior (regime III). A final steady 
state (Regime IV) is reached when all loss processes are 
taken into account. The full sequence of Regimes I - 
IV is most visible when the production rate is not much 
smaller than the interstitial hopping rates. Regimes II 
and III shrink with decreasing production rate, and in the 
limit of very small production rates. Regime I crosses over 



directly into Regime IV. These results were based on pa- 
rameters typical for bcc metals, but the continuum level 
reaction kinetics equally applies to fee metals or other 
crystal structures provided they exhibit similar diffusion 
mechanisms. 

The introduction of vacancy reactions to form immo- 
bile vacancy clusters does not change the general se- 
quence of the scaling regimes, but has a profound effect 
on the population dynamics. In steady state, the density 
of free vacancies and interstitials is reduced relative to the 
non-interacting case, but the total number of vacancies 
nvtot is enhanced. In Regime III, riytot 

and the 

steady state Regime IV is reached at later times. Reason- 
able agreement between rate equations and kMC could 
still be achieved in a point cluster approximation, where 
the capture numbers are size-independent. Of particular 
interest is the size distribution of vacancy clusters, which 
grow and shrink under the diffusive arrival of interstitials 
and vacancies. For irreversible vacancy aggregation, we 
derived a new expression for the cluster size distribu- 
tion, ricirn) — Kl7^™/Km(a\/^?^^)™~"'^. The general form 
of this distribution remains when also allowing vacancy 
detachment from the cluster. 

Immobile interstitial clusters were shown to further in- 
crease the vacancy population. However, interstitial clus- 
ter that are as mobile as single interstitials decay rapidly 
in relevant parameter regimes. We therefore expect that 
nucleation of interstitial clusters through diffusion plays 
a negligible role in the microstructural damage evolution 
of pure bcc metals. It can become important, however, 
if trapping of self- interstitials near impurities occurs. 

Because of their importance in applications and the 
availablitity of detailed molecular dynamics studies, we 
applied the model for parameters suitable for V and a-Fe. 
Within our model, both metals exhibit similar defect ki- 
netics and therefore similar swelling rates at a given tem- 
perature. For production rates and diffusivities in typ- 
ical experimental situations, our calculations revealed a 
sequence of growth regimes for the total vacancy density 
TLytot ~ t^, ~ i^^'^ and ~ The volume swelling rate fol- 
lows the same scaling sequence in the sub-1 dpa regime. 
Interestingly, growth of the vacany cluster density as 
~ i^/^ over several decades in time has been observed 
in the much more detailed kMC simulation of Ref. 
but the origin of this growth law had not been previ- 
ously understood. All our simulations reach the steady 
state regime at damage levels of less than 1 dpa. In the 
experimentally relevant regime of damage levels between 
1 and 100 dpa, it will therefore be more efficient to use 
steady state rate equations rather than explicit kMC to 
predict the long time point defect distribution evolution. 
Nonetheless, it is important to understand the onset of 
damage evolution and the conditions of applicability of 
the steady state theory. 

The fact that experimental swelling rates of V and 
a-Fe are very different is of course an indication that 
our present model does not yet include all relevant phe- 
nomena of radiation damage. One obvious refinement 
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of the models developed here would be a more detailed 
parametrization of the detachment rates (binding ener- 
gies) of interstitial and vacancy cluster. However, the 
intent of the present study has been to focus on gen- 
eral trends rather than quantitative comparisons with 
experiments, and no new physics is expected to ap- 
pear from these additional details. For cascade dam- 
age conditions'^^, however, two other processes not in- 
cluded in this study are known to have a crucial impact 
on the defect evolution (in particular void swelling rates). 
Frenkel pair production is only assumed to be homoge- 
neous for energies right above the displacement thresh- 
old, while higher energies lead to the formation of mo- 
bile interstitial and immobile vacancy clusters during the 
cascade phase. Our present model is therefore most rel- 
evant to low particle energies. Inclusion of intracascade 
clustering processes will change results qualitatively and 
quantitatively, but requires reliable information about 
the cluster size distribution during cascades from MD 



simulations. 

In addition to this "production bias", "absorption 
bias" also usually exists in bcc metals. "Absorption bias" 
leads to a preferred absorption of self-interstitials at sinks 
and is known to be an essential driving force for swelling. 
The origin of this bias, which leads to increased cap- 
ture numbers, are long range elastic interactions between 
point defects and sinks such as dislocations22i2^ or grain 
boundaries'^'*. Extentions of the model to include these 
effects should provide a fruitful topic for future work. 
The combination of kMC and rate theory can then be 
used to determine which parameters are most impor- 
tant in radiation damage, so that atomistic simulation 
resources can be better focused. 
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